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NUMERICAL INTEGRATORS THAT CONTRACT VOLUME 


ROBERT I. MCLACHLAN AND G.R.W. QUISPEL 


Abstract. We study numerical integrators that contract phase space volume even 
when the ODE does so at an arbitrarily small rate. This is done by a split¬ 
ting into two-dimensional contractive systems. We prove a sufficient condition 
for Runge-Kutta methods to have the appropriate contraction property for these 
two-dimensional systems; the midpoint rule is an example. 


1. Introduction. What is a dissipative system? In physics, the term usually refers 
to possession of a scalar function (such as energy) which decreases in time, and one 
speaks of, e.g., the dissipative pendulum, x = — sinx —ex, for which —cosx) = 

—ex^ < 0. (See ^ for some general formulations of such systems.) In dynamical 
systems, it usually refers to a decrease of phase space volume in time, as in the 
“dissipative Henon map” (x, y) h-^ (t/, 1 + 6x — at/^), with Jacobian determinant —b — 
phase space area decreases if |6| < 1. Another example is the famous Lorenz system, 
which contracts volume at a constant rate. In the numerical analysis of ODEs, it has 
been used to describe systems that decrease some norm of the solution, either in the 
sense that < a — h\\x\\^ for some a, 6 > 0, or < 0 for all ||x|| > R > 0 

H- 

In the held of geometric integration, much work has been done in maintaining 
the preservation of a conserved quantity (hrst integral) i, I, 0, the decrease of a 
dissipated quantity (a Lyapunov function) BH, or the preservation of phase space 
volume IB B- Here we look at the missing case, and study how to maintain the 
property of contracting phase space volume. 

Consider the ODE 

(1) X = /(x), X G M” 

with solution x{t) and Jacobian (hrst variation) A{t) = dx{t)/dx{b)) which evolves 
according to 

A = FA, A(0) = /, 

where F{x) = df{x) is the derivative of the vector held /. We have 

— det A = det A tr A j 

= det A tr F 

so that phase space volume contracts, is preserved, or expands when tr F < 0, tr F = 
0, or trF > 0 for all x, respectively. trF is the divergence or trace of the vector 
held /. Strongly contractive systems are those for which there is a 6 such that 

Research at MSRI is supported in part by NSF grant DMS-9701755. 

1 






2 


ROBERT I. MCLACHLAN AND G.R.W. QUISPEL 


tr F < 6 < 0. In this case any consistent numerical integrator will be contractive for 
small enough time step h. Therefore we concentrate on weak contraction (defined 
below), which is a closed property and more difficult to preserve. It turns out that 
requiring contractivity for all h > 0 and all contractive / is prohibitively difficult, 
which leads to the following definition. We consider one-step methods Xn+i = g{xn) 
with Jacobian A = dg{xo). 

Definition 1. The ODE is (weakly) contractive if tr F < 0 for all x. An inte¬ 
grator is (weakly) contractive if for any matrix norm || • || and all L > 0 there is a 
time step h* > 0 such that \ det A| < 1 for all 0 < h < h*, for all x, and for all f 
such that ||F|| < L and trF < 0. 

That is, there might be stiffness problems (for large L, h* might be small), but the 
time step needed to preserve contractivity should not tend to zero as trF —0. Note 
that a contractive integrator as defined here is not necessarily volume-preserving when 
the ODE is, nor is the relative amount of contraction necessarily correct as trF —0. 
These would be true if we added the requirement ln(det A)/htr F —1 uniformly 
as trF —0 uniformly, for all fixed h < h*. The midpoint rule (see Proposition 3, 
below) satisfies this, for example. 

Since there are no known linearly covariant volume-preserving schemes in more 
than two dimensions [0, we expect that the same is true here, and we immediately 
consider systems in two dimensions. 


2. Dissipative schemes in two dimensions. 


Example 2. Euler’s method is not contractive in two dimensions. We have Xn+i = 
Xn + hf{xn) so A = I -\- hF. In two dimensions. 


det A = det 


1 -h hFii hFi2 
hF2i 1 hF22 


1 -\- htiF + hf det F. 


So det d < 1 for all h if det F < 0, and det d < 1 for 


h < 


— tr F 
det F 


if det F > 0, so small contractivity can require a small time step to be captured. 


Note that since det A = det(/ -I- hF) = ])([(1 -|- hAj), where Aj are the eigenvalues of 
F, Euler’s method is contractive in n dimensions on systems with bounded negative 
eigenvalues. We look at this further in Section 3. 


Proposition 3. The midpoint rule, Xn+i = Xn + hf[x), x = {xn + a:n+i)/2, is con¬ 
tractive in two dimensions. 


Proof. We have 


so 


A 


-hF(x 


-1 / 1 

I + -hF(x^ 


det A 


1 -|- he -|- h‘^d 
1 — he K^d 
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where e = |trF(a;), d = jdetF{x). Thus (detA)^ < 1 if 

(l + he + < (l — he + h^d)^ 

or 

e (l + h?c[) < 0 

Since e < 0, this is true for all h if e = 0 (the well-known result that the midpoint rule 
is area-preserving, or symplectic), for all h if d > 0, or for h < 1 /a/^ if d < 0. □ 

Proposition 3 can be generalized as follows. 

Proposition 4. The symplectic Runge-Kutta methods with bi > 0 for all i are con¬ 
tractive in two dimensions. 

Proof. For the terminology, see 0 . Our proof closely follows their proof of symplec- 
ticity. An s-stage Runge-Kutta method is defined by 

S 

( 2 ) Xi = Xn + h'^aijf{Xj), 

i=i 

S 

(3) Xji-\-i Xji h ^ ^ bj f ( Aj) , 

and is symplectic if bibj — biaij — bjaji = 0 for all i and j. Note that in two dimensions, 
A^JA = Jdet A, where J = (_°i J), so we evaluate the left hand side. Let Di = 
df{Xi) = F{Xi)dXi =: FjAj. Differentiating (|]) gives 

A^JA= (^I + hJ2 a) j (/ + h ^ bjDj^ 

i j 

= J+hJ2^i (JDi + Dj J)+h^Yl 
i i,j 

Differentiating (|^) gives 

(4) Aj = I h ^ ^ a^jDj 

j 

or 

JDi = AJ JDi — h cbijDj JDi 
j 

Inserting, 

A^JA = J + h'^^bi (^AjjDi + DfjAfj + A {bibj — biOij — bjOji) DJJDj. 
i i,j 
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The last term is zero because of the assumption on the coefficients bi, aij. Now 
Di = FiAi, so 

A^JA = J + hJ2 biAj (JFi + F'[j)Ai 

i 

= J + h biAj JAi tr Fi 

i 

= j(^l + h bi det Ai tr F^ 

i 

SO 

det A = 1 + h bi det Ai tr Fj. 

i 

From (^), det is bounded and equal to 1 + 0{h). Using fej > 0 and trFj < 0 gives 
the result. □ 

The assumption > 0 is necessary. Suppose there are s = 2 stages with 6 i > 0 
and 62 < 0. Then the vector {bi) lies in the fourth quadrant, and all we know of the 
vector {tr Fi) is that it lies in the third quadrant. In regions where the trace varies 
relatively quickly, the angle between these two vectors can be less than leading to 
{bi) ■ {tr Fi) > 0 and det A > 1 . 

These methods actually preserve area when trF = 0. In fact, this is not necessary 
for contractivity in two dimensions, because we can allow a small amount of “numer¬ 
ical contractivity” even as trF —»• 0 ; away from trF = 0 the inherent contractivity 
of the ODE contributes. It turns out that only methods of order 2, 3, 6 , 7, ... , can 
achieve this. 

Lemma 5. Let R{z) be the linear stability polynomial of a consistent Runge-Kutta 
method. In two dimensions, the method is contractive on linear ODEs if there is a 
M* > 0 such that 

R{u)R{—u) < 1, R{iu)R{—iu) < 1 

for allt) <u < u*. 

Proof. In n dimensions, a Runge-Kutta method on linear problems x = Fx has 
derivative A = R{hF). Therefore det A = = 1 + hti F + 0{h‘^), so the 

method is contractive if trF < 0. If trF = 0 we have to examine detR in more 
detail. In n = 2 dimensions, there are only two such cases: the eigenvalues can be 
{u, —u) or {iu, —iu). This gives the result. □ 

We note that the result also applies to nonlinear problems with 1-stage methods, 
since then F{x) is evaluated at only a single point. 

Proposition 6. Let the method have order p, so that R{z) = + az^^^ + bz^^‘^ -f 

0{zP^^). If 4\{p + 1) and a < 0, or if 4:\{p + 2) and b < a, then the method is 
contractive on linear problems in two dimensions. 

Proof. We expand 

R{z)R{-z) - 1 = e-\azP+^ + bz^^^) + e\a{-z)P+^ + b{-z)P+‘^) P... 

= azP+\l - (-1)^’) + (6 - a)zP+\l + (-If) -F ... 
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The leading term must be negative for z = u and for = m, so it must be a fourth 
power. If p is even, the leading term is so 4|(p + 2) and we need 6 < a; if p is 
odd, the leading term is z^^^ so 4|(p + 1) and we need a < 0. □ 

An example is any 3-stage, 3rd order Runge-Kutta, which has R{z) = 1 + z+^z^ + 

1 „3 

gZ ■ 

This result can be extended to more dimensions. For example, a longer calculation 
shows that Proposition ^ holds with p = 3 in three dimensions. We are not sure how 
it extends to nonlinear systems. It seems that if the eigenvalues of F are varying 
rapidly, contractivity could be lost. 


3. More than two dimensions. For systems in more than two dimensions, we 
generalize the volume-preserving method of Feng and Wang (p[; see also [^). We 
write the ODE as a sum of two-dimensional contractive systems (i.e., ones for which 
Xi = 0 except for two indices i), apply a contractive method to each term, and 
compose the resulting maps with positive time steps. Since contractivity is a semi¬ 
group property, we can build a contractive integrator of order 1 or 2 in this way p[. 
This relies on the following proposition. 

Proposition 7. Any contractive ODE is the sum of two-dimensional con¬ 
tractive ODEs. 


Proof. Consider x = /( t ), F = df. We shall write / in the form /* = Yhj^jLij 
(where dj = d/dxy) 

Let Sij{x) be functions with 

n 

Sij{x) -\- Sji{x) > 0 and Sij{x) = 1 

Li=i 

for all X. Let 

S, = jj 

where any values of the indefinite integrals can be taken. Let 

fi = fi-'^ djSij = /i - ^ j Sij tr F dxi, 


so that 


ti df = f 1 — Sij j tr F = 0. 




Thus, / is traceless and can be written as 




where the matrix A is antisymmetric and as smooth as / ai- Therefore 

( 5 ) = + 


or L = A -I- S'. 
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For an explicit splitting, we take the diagonal elements Sii(x) = 0. Then x = f is 
the sum of the following n(n — l)/2 two-dimensional ODEs: 


x^ dj Lij 
Xj diLji 

Xfc = 0 for /c 7 ^ i, j 


for each pair (f, j) of indices from 1 to n. Each is contractive because each A piece is 
traceless and each S piece has trace [sij + Sji) trF < 0. 

One degree of smoothness is lost in this splitting, because each S piece depends on 
tr F. □ 


An interesting solution is obtained by taking sa = 1/n, Sij = 0 for i 7 ^ j, and 


nLij 



fj dxi A dij 


ti F dxi dxi- 


However, a more practical decomposition is to take the same S but Aij = 0 for 
\i — j\ > 1; this gives the minimum of n — 1 two-dimensional ODEs. 

Although the above proof is constructive, it may be possible to hnd a more con¬ 
venient splitting by ad hoc methods, in some cases leading to an explicit contractive 
integrator. 

Firstly, if / is the sum of integrable contractive vector helds, then their flows can 
be composed to give a contractive integrator for /. For example, the Lorenz system. 


X = 


-a a 0 \ /OO 0 

p —1 0 1 a; -I- I 0 0 —Xi 

0 0 -/3 / \ 0 Xi 0 



1 


is the sum of a linear, contractive part and a Poisson, non-contractive part, each of 
which may be solved exactly, giving an integrator with exactly correct contractivity. 

Secondly, it may be possible to use a simpler method, such as Euler, on some of 
the pieces. Here are some criteria which allow this. 

Proposition 8. Euler’s method is contractive in n dimensions if there is a h such 
that tr(F^) > 6 > 0. This condition is equivalent to ||>S'|p > ||A|p -|- b, where F = 
A + ^, A = -A^, ^ = S^, and || ■ II is the Frobenius (sum of squares) norm. This 
condition is satisfied if all the eigenvalues of F are bounded away from the sectors 
I < 16^1 < in particular, if they are all real and bounded away from zero. 
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Proof. Let A* be the eigenvalues of F. For Euler’s method we have 

In det A = In det(/ + hF) 

= ln]^(l + hXi) 

i 

= ^ ^ ln(l + h\i) 

= h ^2 ~ 2^^ ^2 0{h^) 

i i 

= htiF - ^hHT{F^) + 0{h^). 

If there is a 6 such that tr(F^) > 6 > 0, this is less than 0 for all small enough h, i.e., 
the method is contractive. Splitting F into its symmetric and antisymmetric parts, 

tr(f") = = Y^(S„ + /ly)(S., - A,,) = ||S||2 - ||/lf, 

i,j i,j 

giving the second part of the proposition. Now tr(F^) = if is 

outside the specified sectors, then each real eigenvalue or complex conjugate pair 
of eigenvalues gives a positive contribution to this sum, giving the last part of the 
proposition. □ 

Note that the eigenvalues of elliptic or nearly elliptic fixed points he near the 
imaginary axis—right in the middle of the “bad” sector. Perhaps this was only to be 
expected. 

Experts will recognize the last part of Proposition ^ as the appearance of an order 
star of a Runge-Kutta method (the set {z : |R(2:)| < |e^| where R{z) is the 
method’s linear stability polynomial). For linear problems, or nonlinear problems 
with 1-stage methods, a method is more contractive than the flow of the ODE if hXi 
lies in the order star of the method for each eigenvalue A*. However, this seems rather 
restrictive so we do not explore further. 

Proposition 9. There are explicit contractive integrators. 

Proof. Let / be any contractive vector field with ||F|| < L. Because eigenvalues 
vary continuously and can only become imaginary when two eigenvalues meet, and 
because symmetric matrices have real eigenvalues, there is a symmetric, traceless 
matrix M with distinct eigenvalues such that the derivative of /i := / — Mx has real 
eigenvalues. Let /2 := Mx and split f = fi + f 2 - /i is contractive and admits an 
explicit contractive integrator (e.g. Euler’s method, see Proposition 8); /2 is traceless 
and can be solved explicitly. Composing these maps gives the result. □ 

We close with some open questions we hope to report on in the future. 

1. Are there explicit contractive integrators of any order? (Proposition 9 con¬ 
structs a first order method.) There are if one only demands /inear contractivity. 
The order cannot be increased by composition, because the adjoint of Euler’s 
method—backward Euler—is not contractive for /i. 
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3. 


The present method reduces to the volume-preserving method of Feng and Wang 
0 when the vector held / is traceless. There is another approach to volume¬ 
preserving integration due to Quispel P and to Shang [0 , which does not rely 
on a splitting at all; moreover, it has a generalization to systems preserving 
non-Euclidean measures, which we have not even considered here. Can this 
approach be carried over to the contractive case? 

The splitting used in the proof of Proposition 5 writes f = a+b where tr (da) = 0 
and 6 = 0 when tr((i/) = 0. Are there splittings with the property that b{x) = 0 
when tr((i/(x)) = 0? If so, they could be used for systems in which tY{df{x)) 
changes sign on a compact hypersurface; the interior would then be invariant 
and one could construct an integrator which preserved it and was contractive 
there. This was done for the case of dissipation of scalar functions in 0]. 
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